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Flow of wet granular materials: a numerical study 
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We simulate dense assemblies of frictional spherical grains in steady shear flow under controlled 
normal stress P in the presence of a small amount of an interstitial liquid, which gives rise to capillary 
menisci, assumed isolated (pendular regime), and to attractive forces, which are hysteric: menisci 
form at contact, but do not break until grains are separated by a finite rupture distance. The 
system behavior depends on two dimensionless control parameters: inertial number I and reduced 
pressure P* = aP/('KT), comparing confining forces ~ a^P to meniscus tensile strength Fq = ttFu, 
for grains of diameter a joined by menisci with surface tension P. We pay special attention to 
the quasi-static limit of slow flow and observe systematic, enduring strain localization in some of 
the cohesion-dominated (P* ~ 0.1) systems. Homogeneous steady flows are characterized by the 
dependence of internal friction coefficient /r* and solid fraction $ on I and P*. We also record 
normal stress differences, fairly small but not negligible, and increasing for decreasing P*. The 
system rheology is moderately sensitive to saturation within the pendular regime, but would be 
different in the absence of capillary hysteresis. Capillary forces have a significant effect on the 
macroscopic behavior of the system, up to P* values of several units, especially for longer force 
ranges associated with larger menisci. The concept of effective pressure may be used to predict 
an order of magnitude for the strong increase of p* as P* decreases but such a crude approach is 
unable to account for the complex structural changes induced by capillary cohesion, with a significant 
decrease of $ and different agglomeration states and anisotropic fabric. Likewise, the Mohr-Coulomb 
criterion for pressure-dependent critical states is, at best, an approximation valid within a restricted 
range of pressures, with P* > 1. At small enough P*, large clusters of interacting grains form in 
slow flows, in which liquid bonds survive shear strains of several units. This affects the anisotropies 
associated to different interactions, and the shape of function which departs more slowly from 

its quasistatic limit than in cohesionless systems (possibly explaining the shear banding tendency). 

PACS numbers: 83.80.Fg; 47.57.Qk 
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I. INTRODUCTION 


Over the last decade, constitutive modelling of dense 
granular flows has been proposed [Il,i in terms of in¬ 
ternal friction laws directly applying to normal stress- 
controlled steady shear flows, for which the internal state 
of the material is characterized by a single dimensionless 
number, the inertial parameter I Q. Number I might 
be regarded as a reduced, dimensionless form of shear 
rate 7 = dvi/dx 2 , related to the stress (T 22 normal to 

flow direction x as I = ™ denoting the particle 

mass and a its diameter. The constitutive law relating 
the effective internal friction coefficient, /i*, defined as a 
stress ratio, /r* = oyiIctti-, to inertial number / should 
be supplemented with a similar relation of solid fraction 
$ to / [im. I characterizes dynamical effects, and 
the quasistatic limit is that of vanishing I. In this limit 
of / —>■ 0 , the material is in the so-called critical state 
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of soil mechanics Q , i.e., quasistatic plastic shear flow 
at constant solid fraction under constant stresses and 
effective internal friction /i*. In various experimental and 
numerical studies, the constitutive law, suitably general¬ 
ized, was shown to apply to different grain shapes and 
flow geometries [iS3- A major advantage of the “p*{I) 
and $(/)” approach is its ability to deal with both the 
quasistatic limit and the rigid limit without any diver¬ 
gence or singularity. On regarding inertial number I as 
the sole state parameter in a granular material in shear 
flow, it is implicitly assumed that small contact deflec¬ 
tions due to the finite elastic stiffness of the grains are 
irrelevant - this is the rigid limit. 


In the presence of attractive forces between neighbor¬ 
ing grains, contacts are endowed with a finite tensile 
strength Fq, whence a new dimensionless parameter, P*, 
a reduced pressure comparing the applied confining stress 
P (say, the controlled normal stress value in shear 

u^cr 

flow) to force scale Fq, as P* = - (similarly a “co- 

^0 ,—, 


hesion number” 77 = 1/P* was defined in Ref. Illll). Un¬ 
der small P *, cohesion stabilizes loose structures |l2l. [l3l| , 
which collapse upon increasing P* [l3, [13 • In steady 
shear flow, generalization of rheological laws to the cohe¬ 
sive case involve internal friction coefficient and density, 
functions of both numbers I and P* or 77 [13 ■ 


In wet granular materials cohesion arises from capillary 
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forces due to small liquid bridges ioining particles touch¬ 
ing or in close vicinity to each other [13, [3 ■ The effect of 
such forces has been investigated in quasistatic deforma¬ 
tion (Tgl - f^ , and some of its consequences in terms of mi¬ 
crostructure were discussed [l^l- In the pendular regime 
of saturation US El those bridges are small enough and 
do not merge, so that capillary forces are pairwise addi¬ 
tive. Those attractive forces act as a source of cohesion, 
and are also characterized by a small range and some de¬ 
pendence on intergranular distance, as a liquid meniscus 
might join grains that are not in contact. 

A traditional approach of partially saturated granular 
materials in geomechanics 1231 . w hich has been investi¬ 
gated in recent DEM studies , is to resort to the 

concept of effective stresses, or stresses such that, if ap¬ 
plied to the dry material, would produce the same defor¬ 
mation and plastic flow of the granular skeleton as in the 
ones observed in the wet material. Proposed definitions 
of such an effective stress tensor in the unsaturated case 
generalize the Terzaghi principe [ll applying to satu¬ 
rated media, and involve a correction of the avera ge p res- 
sure related to saturation and capillary pressure |27|. 

On the macroscopic scale, the effect of adhesive forces 
are sometimes described in the quasistatic limit of slow 
flow by the phenomenological Mohr-Coulomb law Bin, 

m, 


ISec. IVIll) . structural anisotropy fSec. IVIIH) . The results 
are discussed and put in perspective in the final, conclu¬ 
sive section 

II. MODEL MATERIAL AND SIMULATION 
SETUP 

We consider a granular assembly composed of N equal¬ 
sized spherical beads of diameter a, made of a material 
with Young modulus E and Poisson ratio v. The con¬ 
tacts are frictional, satisfying Coulomb’s law with fric¬ 
tion coefficient /i. The granular flow is set by imposing 
a uniform shear rate 7 to a rectangular parallelipipedic 
cell with edge lengths (La)i<a< 3 - In order to avoid wall 
effects periodic boundary conditions are used in all three 
directions. The periodicity, in the direction of the flow 
gradient (direction 2), is applied with the Lees-Edwards 
procedure [13] and in the two other directions the bound¬ 
ary condition is simple periodic. The system size L 2 is 
allowed to fluctuate in order to keep the normal stress 
(T 22 constant equal to a prescribed value P while Li and 
L 2 are fixed Q. 

A. Force model 


= c + Mto-22, (1) 

characterized by macroscopic cohesion c and internal fric¬ 
tion coefficient p,\. 

The present paper investigates the constitutive laws 
applying to wet model granular materials, in the pendu¬ 
lar regime, and discusses the influence of capillary effects 
on macroscopic behavior and microstructure. Similarly 
to refs. [l 3 , [l 3 |, the rheology and micromechanical as¬ 
pects are studied for varying P* and / values (with spe¬ 
cial emphasis on the quasistatic limit of / —>■ 0). As in 
dry granular systems and in previous studies on 2D cohe¬ 
sive materials the material rheology is described in terms 
of apparent friction coefficient (stress ratio) p* and solid 
fraction d) as functions of I and P *, and the applicability 
of a Mohr-Coulomb relation is tested. Rheological and 
microstructure features as normal stress differences and 
formation of large clusters bonded by liquid bridges are 
also investigated. 

In the following, we first introduce (Sec. |IT|) the mi¬ 
croscopic ingredients of the model material, and report 
then, in Sec. mil on the conditions in which homogeneous 
steady states are observed in shear flows, enabling ma¬ 
terial constitutive laws to be deduced. Such laws are 
measured, depending on the relevant dimensionless pa¬ 
rameters and on some features of the microscopic model, 
in Sec. El Next, in Sec. El we investigate the role of 
capillary forces and distant interactions in the material 
rheology, and revisit the traditional concepts of effective 
stress and Mohr-Coulomb cohesion. Additional stud¬ 
ies of microstructural and micromechanical aspects fol¬ 
low: force distributions (Sec. I VII) . agglomeration effects 


Elastic and frictional forces are jointly implemented 
in contacts as in Ref. [3l|, in which a simplified Hertz- 
Mindlin-Deresiewicz force model is used for the elasto- 
plastic contact behaviour. This model combines the nor¬ 
mal Hertz force F^r, depending on contact deflection h 
as 


Fn = 



( 2 ) 


in which we introduced notation E = E/(I — with a 
tangential elastic force Ft-, to be evaluated incrementally 
in each time step of the simulation. The simplification of 
tangential elasticity adopted is that of Ref. [^ , involving 
a constant ratio {2—2v)/{2 — v) of tangential {Kt) to nor¬ 
mal {Kn) stiffnesses in contacts, both depending on Fv, 

dFjsf 

as, from @, one has Km = ,, oc F^/^. Caution should 

dh 


be exercised to avoid spurious creation of elastic energy 
with such laws, and Kt should be suitably rescaled in 
cases of decreasing normal force and deflection. For the 
details of the elastic model, for the enforcement of the 
Coulomb condition IIFtH < and for the objective 

implementation of the force law, with due account of all 
possible motions of a pair of contacting grains, the reader 
is referred to [M |. 

An estimate of the typical contact deflection under con¬ 
fining stress P defines a dimensionless parameter, stiff¬ 
ness number k [^ . such that h/a oc k~^. For a Hertzian 
contact, one may use [3l| 




( 3 ) 
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Two values of k, 8400 and 39000, used in this study, 
respectively correspond to glass beads with E = 70 GPa 
and 1 / = 0.3 under pressures 100 kPa and 10 kPa, as 
in Ref. Finally, the force model of [3l| which is used 
here may also comprise a viscous damping term opposing 
normal relative motion of contacting grains, chosen to 
correspond to a restitution coefficient close to zero in 
normal collisions. We do not comment this feature, as it 
was shown 0,0 to have very little influence in the slow 
shear flows of the present study. 

The presence of a small amount of an interstitial wet¬ 
ting liquid introduces additional capillary forces, trans¬ 
mitted between contacting or neighboring grains by a 
liquid bridge, or meniscus, as sketched in Fig. [H We 



FIG. 1: (Color online) A meniscus between two spherical 
grains of diameter a = 2R, with distance h between solid 
surfaces, filling angle i/J, contact angle 9. 


consider a perfectly wetting liquid, with contact angle 
9 equal to zero. As in Ref. [l^, we assume that the 
menisci only form for touching particles, but breaks for 
gaps larger than a certain rupture distance Dq, as ob¬ 
served in Dq relates to meniscus volume V as 

Do ~ 

For the attractive force between particles separated 
by distance h < Dq we adopt the Maugis approxima¬ 
tion [ 3 ^, which is appropriate for small enough menis¬ 
cus volume, for its simplicity. The maximum attractive 
force (tensile strength) is reached for contacting particles, 
and equal, according to this model, to Fq = ttoF (F is 
the liquid surface tension) independently of the meniscus 
volume. The capillary force varies with the distance h 
between particle surfaces as 

{ -Fo h<0 

1 ^ ] 0 < ^ < -Do (4) 

0 > Dq 

h < 0 corresponds to an elastic deflection of the particles 
in contact. This formula is a simpler, analytical form of 
the toroidal approximation with the “gorge method” 
for the capillary force in a meniscus, which describes the 
meniscus as limited by circular arcs in a plane containing 
the two sphere centers. An alternative form was given by 
Willett et al. [s^ , while Soulie et al. [ 2 ^ HI, proposed 
a parametrized numerical solution. Fig.[2]compares func¬ 
tions F'^’^P{h) according to Maugis and to Soulie et al.. 
Some comparisons between several formulae and experi¬ 
ments are given in [s^. 



FIG. 2: (Color online) Force law for two different 

meniscus volumes, according to the Maugis model and to the 
Soulie formula. 


B. Saturation range of pendular state 


The morphology of partially saturate d g ranular mate¬ 
rials depends on the liquid content [H, El- The present 
study, like a number of previous ones [13,1^ , is re¬ 

stricted to the pendular state of low saturations, in which 
the wetting liquid is confined in bonds or menisci join¬ 
ing contacting grains. Liquid saturation S is defined 
as the ratio of liquid volume fl/ to interstitial volume 
D,y. It is related to meniscus volume F, solid fraction 
<i> = 1 — Fly/D. and wet coordination number z (the aver¬ 
age number of liquid bonds on one grain) as: 


fli _3z 4> V 

fly TT 1 — <!> 


(5) 


In our study, we fix the value of meniscus volume V. 
Such a choice does not conserve the total liquid volume, 
which is proportional to the varying coordination number 
2 of liquid bonds. Its consequences have to be assessed, 
and we shall check that the results are not significantly 
affected, within the range of investigated material states. 

The pendular state to which our model applies is only 
relevant in some limited saturation range. On the one 
hand, a minimum liquid volume is necessary for menisci 
to form at contacts, as the liquid will first cover the grain 
surface asperities. This minimum saturation 5'niin for 
bridges to form might be roughly estimated upon intro¬ 
ducing a roughness scale S, assuming a layer of thickness 
S covers the surface of the grains, as 


(1 - $)a' 


( 6 ) 


For $ = 0.5 and S ~ 10“^a the minimum value for sat¬ 
uration is of the order of 10“^, as observed in experi¬ 
ments 03- Using and typical values of z (5 or 6) 
and $ (say, 0.6), this sets a lower bound to meniscus 
volume, of order 10“"^a^. On the other hand, the upper 
saturation limit for the pendular state corresponds to 
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the merging of the menisci pertaining to the same grain, 
which, considering a triangle of spherical grains in mutual 
contact, happens as soon as the filling angle (se e Fig. [T]) 
reaches tt/G. The analytical formula for V |34l |. within 
the toroidal approximation, as a function of (p (setting 

h = 0, and 0 = 0 ), then yields — 8 . 10 “^, correspond- 

ing, using ®, to a maximum saturation between 0.05 
and 0 . 1 , similar to experimental observations [ 13 , [ 11 . 


C. Choice of parameters 

Tab.U gives the values of parameters employed in our 
simulations. While stiffness number k and friction coef¬ 
ficient p, are kept fixed, reduced pressure varies from the 
dry case P* = oo down to the lowest value 0.1, for which 
cohesive effects are strong, while the investigated range 
of I values allows us to approach the quasistatic limit 
with some accuracy, as well as assess the effects of iner¬ 
tia in faster flows (although rapid, strongly agitated flow 
are not studied here). The meniscus volume is chosen as 
V = 10“^a^ in most simidations. A few tests are car¬ 
ried out with different values (as given within brackets) 
of the number of particles, the stiffness number and the 
preset meniscus volume. Taking F = 7.3 x 10“^ J.m“^ 



8400 (occasionally 39000) 

0 

0.3 

N 

4000 (8000) 

I 

from 10”"^ to 0.562 by factors of VTO 

P* 

0.1 ; 0.436 ; 1 ; 2 ; 5 ; 10 ; 00 

V/a^ 

10“® (10"^ ; 5 X 10"® ; 2 x 10"'‘ ; 10"®) 


TABLE I: List of parameter values: N particles of diameter a, 
interacting with friction coefficient p, forming menisci of vol¬ 
ume V at contacts, are subjected to normal stress-controlled 
shear flow for which inertial nnmber I, reduced pressure P* 
(evaluated with normal stress a^^) and stiffness parameter k 
take values as prescribed. Attractive forces fall to zero at 
distance Do = . 

for water, and a = 0.1 mm, P* = 1 corresponds to 
CT 22 = (ttF/o) ~ 2.3 kPa - the pressure, under gravity, 
below a granular layer with a thickness of a few tens of 
centimeters. 


as averages over time series. Stresses are measured using 
the standard formula, for all coordinate index pairs a, /3 



i i<j 


(7) 


involving a kinetic contribution with a sum over grains 
f, of velocities Vi and a sum over all pairs with center-to 
center vector r^, interacting with force Fij, denoting 
the sample volume. 

The evolution of solid fraction $ with strain 7 is shown 
in Fig. [3] $ decreases until it approaches its steady state 
value for 7 > 5 in this case. Fig. |4] shows the evolutions 
of (7^2 ®’i 2 with 7 . (Note that is negative with our 

sign convention). We thus check that normal stress ^22 is 
well controlled since its fluctuations about its prescribed 
value P are very small. Shear stress exhibits a fast 
increase and and overshoot at small strain and then de¬ 
creases, approaching its steady state value, after a few 
strain units, over a strain interval (a few units) similar 
to the one corresponding to the transient evolution of $. 
Stresses and solid fraction fluctuate in the steady state, 
and a careful evaluation of measurement errors on their 
time averages is required (especially for shear stress, for 
which fluctuations levels reaching about 20 % of the mean 
value are apparent in the example of Fig. n. We use the 
blocking technique of Ref. [1^ to estimate error bars on 
averages over finite time series. 



FIG. 3: (Color online) Solid fraction versus shear strain 7 . 
Time series is obtained with P* — 1, I = 0.1 and N = 4000 
when the rupture distance is Do = 0 . 1 . 


III. HOMOGENEITY AND STATIONARITY 
A. Steady states, macroscopic measurements 

Starting from a dense initial configuration, with solid 
fraction close to the random close packing value (‘f’ucp — 
0.64), we impose a constant shear rate 7 and wait until a 
steady state is reached before measuring constitutive re¬ 
lations for stresses and solid fraction, which are identified 


B. Shear localization 

1. Velocity profile 

Instantaneous velocity profiles u®(= {v^{x 2 ))s) are 
computed on averaging particle velocities along the mean 
flow direction within slices of thickness 0 . 01 L 2 . We ob¬ 
serve a strong localization of the flow for the smallest 
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FIG. 4: (Color online) Shear stress \o-i 2 \ (lower curve, red, left 
axis) and normal stress cTjj (upper curve, blue, right axis) 
versus shear strain 7 . Note the different scales on left and 
right axes. Time series obtained with P* = 1, I = 0.1 and 
N = 4000 when the rupture distance is Do = 0.1. 

value of the reduced pressure, P* = 0.1, for both slow 
and flows. As represented in Fig. [5^, the velocity gra¬ 
dient, initially uniform, gradually concentrates within a 
shear band of thickness H < 3a which may move verti¬ 
cally but persists for all values of strain 7 > 250. Local¬ 
ization tendencies in slow flow of dry granular materials 
are sometimes reported flain, although, for uniform 
strain rates, usually not observed as an enduring, sys¬ 
tematic phenomenon. In the present study persistent lo¬ 
calization profiles are also detected in nearly quasistatic 
flows, with a shear band thickness between 5a to 10a. 
However, for the intermediate values of the inertial num¬ 
ber ( 10 “^ ^ ^ ^ 10 “^) this effect diminishes and the 
strongly localized profiles are less frequent. 

For all P* > 0.436, localization is not frequently ob¬ 
served and temporary, even in the quasistatic limit. The 
velocity profiles for P* = 0.436 and I = 10“^, as repre¬ 
sented in Fig. [Sj), are nearly linear and on average the 
flow is homogeneous. 


(b) 



FIG. 5: (Color online) Velocity profile for P* = 0.1, I = 0.178 
(a) and for P* = 0.436 and 7 = 10~® (b) at different shear 
strain values. 

instance, when I = 10“^, 4)® decreases from 0.47 to about 
0.4 inside the sheared zone of thickness H k. 7a. 


3. Deviation from linear profile 


2. Local solid fraction 

Similarly we can calculate a solid fraction profile, on 
averaging the solid contents of slices orthogonal to the 
velocity gradient (splitting the volume of one grain be¬ 
tween different slices if necessary). Fig. [5] shows the ve¬ 
locity (u®) and solid fraction ($®) profiles for two differ¬ 
ent values of shear strain, 7=1 and 7 = 352, which 
belong to the simulation of Fig. [S^.. We see that in the 
homogeneous flow the distribution of mass in the sys¬ 
tem is almost uniform, but when the localization occurs 
4’® strongly decreases within the shear bands, to a value 
below 0.2. It slightly increases outside the shear band, es¬ 
pecially in its vicinity. Slighter decrease of density within 
thicker shear bands in quasistatic flow is observed. For 


The deviation from the linear profile is characterized 
by parameter A: 


L ,/2 


/ ivAx2)-iX2fdX2. 


( 8 ) 


- L ,/2 


The normalization by —- ensures a maximum value 
12 


A = 1 in the case of a perfect localization within a plane, 
as if two solid blocks were sliding on each other. As de¬ 
fined in Eq. ®, parameter A is not affected by a vertical 
shift in the velocity profile, due to the periodic boundary 
condition in direction X 2 . If the strain rate is homoge¬ 
neous within a shear band of thickness H and vanishes 
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FIG. 6 : (Color online) Velocity profiles (lower axis in red) and 
local densities of grains (upper axis in blue) for two configu¬ 
ration samples with shear strain 7 = 20 (a) and 7 = 352 (b) 
when P* = 0.1 and I = 0.178. Average solid fraction ($'’), 
is shown with a vertical dashed line in blue, with a value of 
0.44 in (a) and 0.49 in (b). 


outside, one observes 

A = (1-^)2 (9) 

Fig.[7]is the plot of A as a function of strain 7 correspond¬ 
ing to the same simulation as in Fig.[^. It initially shows 
small fluctuations near zero, and suddenly increases near 
7 = 250 when the velocity gradient localizes in a shear 
band. 


4- Occurence of shear banding 

Large values A > 0.8 for P* = 0.1 in the faster flows 
(/ > 0.17) indicate strong localization in this range. At 
/ = 0.1 A drops down to small values, typically below 
0 . 1 , but in the quasistatic limit it increases again: at 



FIG. 7: (Color online) Deviation from linear profile, A( 7 ), 
versus strain 7 for P* — 0.1, I = 0.178. 


I = 10“^ it mainly fluctuates between 0.4 and 0.8. For 
P* > 0.436 the shear rate is much more homogeneous. 
A almost vanishes in faster flows, increases somewhat in 
the quasistatic limit, but rarely exceeds 0 . 2 , even for the 
smallest inertial numbers. 

Simulations carried out with a larger stiffness number 
{k = 39000), for the two smallest values of P* and for all 
values of I in Tab. HI do not record any significant influ¬ 
ence of K on the homogeneity of the flow. The influence 
of the sample size is studied by simulating some samples 
with height twice as large as in the standard sample, 
containing 8000 grains, with P* = 0.436 and different 
values of I. The size dependence in formula |9] implies 
then larger values of A, should the shear strain tend to 
localize, temporarily or permanently, within a region of 
fixed thickness. In our tall, 8000 grain systems, as the 
quasistatic limit is approached, A reaches peak values 
above 0.4 but continuously evolves and no persistent lo¬ 
calization pattern is detected. 

Consequently, our results reveal a strong localization 
tendency at P* = 0.1 for both small (below 0.03) and 
large (above 0.3) values of the inertial number. We per¬ 
formed some measurements at P* =0.1 for intermediate 
values of /, over strain intervals for which values of in¬ 
homogeneity parameter A averaged below 0.1, as in the 
first part of the graph of FiglT] We limited our system¬ 
atic studies to homogeneous flows, for larger values of 
P *, for which no evidence of enduring localization effects 
is observed. 


A systematic fluid depletion in shear bands was re¬ 
ported in 4^ - this requires a model for liquid migra¬ 
tion between menisci, which we did not introduce in the 
present study. We limit our results to the issue of whether 
shear banding occurs, and focus on homogeneous flows. 
Shear banding and inhomogeneous flows certainly de¬ 
serve more detailed studies, which we plan to pursue in 
another publication. 
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IV. MACROSCOPIC BEHAVIOUR AND 
CONSTITUTIVE RELATIONS 

Focussing on homogeneous steady states (excluding 
too small P* values, from Sec. IIII Bl) we now deduce 
macroscopic constitutive relations from the simulations. 


A. Shear stress and solid fraction 


Friction coefficient /r* and solid fraction depending 
on / for various P* values, are shown in Fig. [5] for the 
parameter choice adopted in most simulations. We fit the 


(a) 



I 

(b) 



I 

FIG. 8: (Color online) Macroscopic friction coefficient fi* (a) 
and solid fraction <1? (b) versus inertial number I for different 
values of reduced pressure P* {Do = 0.1a = 


following power law functions to those data, denoting as 
and $0 the quasistatic limits (critical state) values 
of the macroscopic friction coefficient and of the solid 
fraction: 


= mS + c/“ 

$-1 = $-1 + e /'^ 


( 10 ) 


Tabs. El and uni give the values of the fitting parameters 
introduced in Eqs. (TU] While the increase of and the 
decrease of $ as functions of / are familiar trends, similar 
to observations made with dry grains [I|, , some other 

features are remarkable. The quasistatic limit is quite 
nearly approached for I < 0.01, and is strongly influ¬ 
enced by capillary forces. Internal friction coefficient p*, 
compared to the dry, cohesionless value (0.335±0.001), 
already shows a notable increase at P* = 10, reaching 
values as high as 0.6 for P* = 1 (i.e., as cohesive and 
confining forces are of the same order), and nearly 0.9 
for P* = 0.436, about 2.3 times the cohesionless value. 
Our partial results for P* = 0.1, measured in reason¬ 
ably homogeneous flows, indicate ~ 1.6 for I = 10~^. 
Meanwhile, the material becomes looser, with reach¬ 
ing values that cannot be observed without cohesion in 
quasistatic conditions. 


p* 

Po 

a 

c 

0.436 

0.867 ±0.003 

0.76 ± 0.04 

0.36 ±0.01 

1 

0.607 ± 0.003 

0.79 ±0.05 

0.41 ±0.02 

2 

0.477 ± 0.005 

0.76 ± 0.05 

0.45 ± 0.02 

5 

0.391 ± 0.006 

0.73 ± 0.05 

0.49 ± 0.02 

10 

0.367 ± 0.004 

0.79 ± 0.04 

0.51 ±0.02 

CX2 

0.335 ± 0.001 

0.84 ± 0.02 

0.62 ±0.02 


TABLE II: Parameters of the fit of function p* (/) by Eq. 1101 
for different values of P*. 


p* 

$0 

V 

e 

0.4360 

0.5243 ± 2.lO""' 

1.73 ±0.05 

0.497 ±0.017 

1 

0.5559 ± lO""' 

1.34 ±0.012 

0.512 ±0.005 

2 

0.5726 ± lO""' 

1.21 ±0.01 

0.547 ±0.003 

5 

0.5851 ± 10"^ 

1.12 ±0.01 

0.580 ± 0.003 

10 

0.5900 ± 10"^ 

1.09 ±0.01 

0.594 ± 0.004 

OO 

0.5970 ± 10"'‘ 

0.96 ±0.015 

0.562 ± 0.008 


TABLE III: Parameters of the fit of function 4>(/) by Eg. 1101 
for different values of P*. 


Such a strong influence of cohesive (capillary) forces 
contrasts with the results of Refs. [iTI. [iq , in which simi¬ 
lar deviations between cohesionless and cohesive systems 
are not observed until P* decreases to much lower val¬ 
ues, of order 0.01. Such 2D results were however obtained 
with a different attractive force law, of vanishing range 
beyond contact. The origins of the strong rheological ef¬ 
fects of capillary forces are investigated in the following. 



























B. Normal stress differences 


n on 


(a) 


The first and the second normal stress differences are 
defined as 



( 11 ) 


Note that those definitions coincide with the one used in 
complex fluid or suspension rheology , but that we use 
the opposite sign convention for normal stresses. Signs 
of Ni and N 2 should thus be reversed for comparisons to 
this literature. 

Most often, considering dense flows of dry granular ma¬ 
terials, those differences, deemed small, are ignored or 
neglected [ 2 ^. We find it worthwhile to record their val¬ 
ues nevertheless, since, as shown in Fig. [HI where Ni and 
N 2 are plotted versus I for different values of P*, they 
are strongly influenced by capillary forces. The first nor¬ 
mal stress difference is very small in the quasistatic limit 
and for large values of the reduced pressure. It increases 
with I and for decreasing values of P*, going through 
a transition from small negative values to positive val¬ 
ues near I = 0.03, for P* > 2. 7Vi variations with / are 
nearly parallel for different P* values. The second normal 
stress difference N 2 also increases for faster flows and for 
decreasing reduced pressure P*. In the quasistatic limit, 
it is considerably larger then A^i. 


C. Sensitivity to capillary force model and 
saturation 



1. Capillary force model 

We tested the effect of the capillary force model by re¬ 
placing the Maugis approximation, Eq. |4l with the more 
accurate parametrized capillary force law proposed by 
Soulie et al. [10, for V = Although the dif¬ 

ference in the force models is appreciable on a plot of 
^Cap ygj-g^g (with the Soulie force about 10% smaller 
at contact, see Fig. ED, the macroscopic results are very 
close: the difference in stress ratio fi* and solid fraction 
$ increases with I but does not exceed 2%. 


2. Meniscus volume and force range 

Changing the meniscus volume amounts to changing 
the distance at which the attractive force vanishes, rup¬ 
ture distance Dq = as well as the gap dependence 
of the capillary force (Fig. [2|). Fig. [TUI shows in¬ 

ternal friction coefficient pi* {I) to be strongly sensitive to 
meniscus volume for the lowest P* values. For a menis¬ 
cus volume of 10“®a^, as compared to the standard value 
10“^a^, fi* decreases by about 20%. Actually, for such a 
small meniscus volume, the decay of the attractive force 
(Fig. [U|) is so fast that, as we checked, results are hardly 


FIG. 9: (Color online) First (a) and second (b) normal stress 
differences as functions of 1 for different P* values. The same 
symbols and colour codes apply to both figures. 


changed on setting the meniscus rupture distance to zero. 
The effect on the solid fraction <I> remains small. To ex¬ 
plore the rheological properties throughout the pendular 
regime, we varied the meniscus volume, and recorded the 
solid fraction and the friction coefficient in the quasistatic 
limit for the smallest studied P* value, as indicated in 
Tab. lIVi thus fully covering the corresponding satura¬ 
tion range (see Sec. IIIBI) . Saturation S, by relation ([5|), 
is related to the wet coordination number, z, whose val¬ 
ues are also provided in the table. While the change 
in solid fraction does not exceed 0.01, the variation of 
the macroscopic friction coefficient is about 20% in the 
pendular regime (up to 50% upon extending the numer¬ 
ical study to unrealistically small menisci, V = 10“®a^). 
We therefore predict a moderate variation of rheological 
properties within the simulated pendular regime of the 
partially saturated granular assembly. Returning to the 
basic assumptions of our model, one of its drawbacks is 
that it ignores liquid volume conservation. Within the 
granular sample, the total liquid volume is proportional 
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FIG. 10: (Color online) Macroscopic friction coefficient /r* (a) 
and solid fraction >I> (b) versus inertial number I for different 
values of P* and meniscus volume V (with Do = 


Vja^ 


Z 

-F 

d* 

10 "^ 

7.137 X 10"^ 

6.863 

0.520 

1.071 

5 X 10"® 

3.418 X 10"^ 

6.556 

0.522 

1.003 

10 "® 

6.305 X 10"® 

5.970 

0.524 

0.875 

2 X 10"'‘ 

1.075 X 10"® 

5.534 

0.525 

0.787 

10 "® 

5.539 X 10"® 

4.836 

0.530 

0.661 


TABLE IV: (Color online) Effect of meniscus volume or 
saturation level on different parameters for I = 10“^ and 
P* = 0.436. 


to coordination number z. As z varies with I, we should 
in principle correct the meniscus volume to maintain a 
constant product zV for different shear rates. However, 
the V dependence of macroscopic properties is so slow 
{fi* varies by 20% as V is multiplied by 20) that the re¬ 
sulting correction on V (as z changes, typically, from 6 
to 4 at most) should not exceed 2%. 


3. Hydraulic hysteresis 

Another feature of meniscus model the role of which 
should be explored is the hysteresis of the attractive force, 
which appears at contact, and vanishes at distance Dq. 
Given the results of Tab. lYl one may expect a strongly 
enhanced influence of distant interactions (as reported, 
e.g., in Ref. [l^). Fig. [TTl compares internal friction and 
solid fraction, for different values of P* and I, in the stan¬ 
dard, hysteretic model, and without the capillary force 
hysteresis, assuming menisci to appear as soon as non¬ 
contacting grains approach below distance Dq. With the 


(a) 



(b) 



FIG. 11: (Color online) Macroscopic friction coefficient fi* 
(a) and solid fraction 4> (b) versus inertial number I for P* = 
0.436 and Do — 0.1, with and without capillary hysteresis. 

new rule of meniscus formation, strongly increases, es¬ 
pecially for small values of I. The internal friction yL *, for 
I ~ 0.1, is close to the standard case, but larger values 
are obtained as I decreases. Even for the smallest values 
of / investigated (/ = 0.001), the material properties still 
depend on shear rate and no proper critical state appears 
to be approached in our simulations. The decrease of /i* 
as a function of I in interval 0.001 < / < 0.01 should trig- 
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ger shear-banding instabilities, as discussed in [11,14^. A 
slightly decreasing trend of versus I was also appar¬ 
ent in FigfTOl for very small Dq. For the standard value 
V = 10“^a^ adopted in this study (as one for which lab¬ 
oratory observations should be possible), the friction co¬ 
efficient does increase with I, albeit slower and slower as 
P* decreases (see coefficient c in TabHIl). The stabilizing 
effect of this growing variation is weaker as cohesion gets 
stronger, which is consistent with the systematic shear 
banding behaviour at P* = 0.1, and might be jeopar¬ 
dized on tampering with the capillary force model. 


V. RHEOLOGICAL EFFECT OF CAPILLARY 
FORCES 


We now seek to explain the strong influence of capillary 
forces on the macroscopic material rheology. The roles 
of different interactions, in the force network and in the 
stresses are investigated. We first collect information on 
coordination numbers and neighbor distances ISec IV Al) . 
Simple relations to average forces are recalled in Sec. lVBl 
We split the stresses into several contributions, in order 
to appreciate the importance of different types of forces. 
This decomposition (Sec. IV Cl) suggests an attempt to 
relate the rheology of wet grains to that of dry ones, 
in terms of some ’’effective pressure” approach, in the 
quasistatic limit, which we present in Sec. IV PI 


A. Coordination numbers and near neighbor 
distances 

Fig.[T2]shows the variations, with inertial number /, for 
different P* values, of coordination numbers Zc, for pairs 
of grains in contact, and Zd, for pairs of grains attracting 
each other without contact at a distance lower than Dq. 
The average number of contacts per grain, Zc, decreases 
for larger inertial numbers, as previously observed in co¬ 
hesionless systems [H, 0 and in cohesive ones [l^ , slower 
for smaller P*, as in [l^ too. Zc also increases as P* de¬ 
creases at constant I, as previously observed as well [l^ . 
Note that this latter trend is opposite to that of the solid 
fraction (Fig. [5]). As the importance of adhesion, relative 
to confinement stresses, increases, looser systems are ob¬ 
tained, yet better coordinated. Grains tend to stick to 
one another, and may form loose aggregates, as in static 
or quasistatically compressed assemblies, for which little 
correlation is also observed Bin between density and 
coordination number. On the other hand, the variations 
of the coordination number of distant interactions, z^, 
with both parameters I and P *, are in the opposite direc¬ 
tion to those of Zc . As / increases, so does Zd- contacting 
pairs tend to separate, but some remain bonded by liq¬ 
uid bridges. And for stronger cohesion (smaller P*), Zd is 
correlated with the system density. The faster approach 
to quasistatic limit at smaller values of P* is apparent in 
both figures. The fraction of rattlers (beads carrying no 


force i) in non-cohesive systems is about 5%. In the co¬ 
hesive case, due to the attractive forces, nearly all of the 
particles are bonded to others and the number of rattlers 
tends to zero, as observed in 2D simulation of cohesive 
powders B; Bl • tends to compensate the changes of 


(a) 



I 

(b) 



I 


FIG. 12: (Color online) Coordination numbers: (a) of con¬ 
tacts, Zc', (b) of distant interactions, Zd- 

Zc, SO that the total coordination number z = Zc -b Zd, 
throughout the investigated range of / and P* values, 
exhibits rather small variations (see Fig. [T^ . Within the 
investigated parameter range, the maximum change in z, 
between 6.8 and 4.8, corresponds to a correction of inter¬ 
nal friction fi *, should we change the meniscus volume to 
maintain the total liquid volume constant, of about 5% 
(see Tab. IlYl). The contact coordination number does not 
change much with the force range or the meniscus vol¬ 
ume. Setting Dq = 0 (instead of the standard value 0.1 
used in the present study) or decreasing the volume of 
the meniscus from its standard value V = down 

to merely leads to a small decrease of Zc, from 5 

to 4.7 in the quasistatic limit, when P* = 0.436. How¬ 
ever, it has a strong influence on Zd. Compared to the 
standard case, for P* = 0.436 and small values of I, it 
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FIG. 13: (Color online) Total coordination number z = Zc + 
Zd- 


decreases from 0.9 down to a value below 0.3 when we set 
Dq = 0.01, or down to about zero when we set V = 10“®. 

It is interesting to compare the number of distant, in¬ 
teracting pairs to the total number of neighbor pairs at 
distance below Dq. The coordination number, z(/i), of 
neighbor grains at distance below h (such that 2(0) = Zc) 
grows with h as depicted in Fig. [TH corresponding to 
/ = 10“^ (quasistatic limit). z{h), like the contact coor¬ 
dination number, is a decreasing function of P* for small 
h/a (below about 2.5 x 10“^, see the insert in Fig. [T4l) . 
It increases with P*, like the density, beyond that dis¬ 
tance. In denser systems grains have more neighbors on 
average, but this is only true if neighbors at some dis¬ 
tance are included in the count, and does not apply to 
contacts (a situation reminiscent of some observations in 
static packings of cohesionless grains [3l|). Up to menis¬ 
cus rupture distance Dq, equal to 0.1a in the present 
case, each grain has on average z{Dq) — Zc non-contacting 
neighbors, among which zq are joined by a liquid bridge. 
Values of ratio zq/{z{Do) — Zc) for different P* and / are 
given in Tab. [V] The proportion of the neighbors within 
range Dq that are bonded by a liquid bridge varies from 
0.61 to 0.71 for / = 10“^ and between 0.68 and 0.79 for 
I = 10“^ - slightly larger than the proportion ^ 50% 
reported by Kohonen et al. in static grain packs. 



I = 

10 “® 

7 = 10"^ 

p* 

Zd 

Zd/{z{Do) - Zc) 

Zd 

Zd/(z(I)o) - Zc) 

0.436 

0.923 

0.609 

1.28 

0.681 

1 

1.27 

0.650 

1.80 

0.723 

2 

1.54 

0.675 

2.24 

0.751 

5 

1.83 

0.701 

2.71 

0.776 

10 

1.97 

0.712 

2.93 

0.787 


TABLE V: Distant coordination number Zd and proportion of 
pairs within distance Dq joined by a meniscus, Zd/{z{DQ) — 
Zc), versus P* for two different values of I. 



FIG. 14: (Color online) Goordination number of neighbor 
grains versus interparticle distance h, for different P* values 
and for I — 10“®. Inset shows detail at small h. 

If the meniscus forms as soon as grains approach at dis¬ 
tance Dq, rather than at contact, the number of contacts 
hardly changes {zc increases by about 5% for P* = 0.436 
in the quasistatic limit), but the increase in the num¬ 
ber of menisci is larger than expected from the data of 
Tab.ig from a simple count of pairs within range Dq : Zd 
is multiplied by 1.7 at small P* and I. 

B. Pressure and average normal forces 

From Eq. (3 neglecting the deflection of contacts in 
comparison to grain diameter a, and ignoring the ki¬ 
netic term, one may relate [3l| the average pressure, 
V = trCT/3, to the average normal force {F^) for all 
interactions, and to the average, {F^h)d, over pairs in 
distant interactions, of the product of force by distance 
h < Dq: 

V = ( 12 ) 

7ra TTO'^ 

. V . 

Due to normal stress differences, ratio — is only slightly 

^ 22 

different from 1 (about 0.95) at small /. We checked that 
formula (TT^ is very accurate for all P* values, and found 
its second term to be negligible, contributing less than 2% 
of the pressure. 

C. Contributions to stresses 

The contribution of the kinetic term in Eq. [7]to stresses 
is quite small. Even for the fastest flow in our simulation 
(/ = 0.562), this contribution does not exceed 2% of the 
shear stress or 5% of the normal stress components, and 
for 1 = 0.178 it is nearly zero for all stress components. 
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Therefore, in this section we only discuss the contribu¬ 
tions of forces to the stress components, for the different 
values of the control parameters, P* and I. These contri¬ 
butions may be split in different ways, on distinguishing 
different forces. 


1. Contact forces and distant capillary attraction 

First one may consider the total stress as a sum of the 
contributions of the contacts and of the distant interact¬ 
ing pairs, as 

CTa/3 =<;3 + '^a/3- (13) 

Our results show that the contribution of contact forces 
dominate in the shear stress. It is larger than 90%, re¬ 
gardless of the values of P* and I. The contribution of 
distant interactions to , as represented in Fig. [15] a, 
is not negligible and increases for smaller values of P*. 
However, it hardly exceeds 10% of the total shear stress. 

The contribution of distant interactions to is dis¬ 
played in Fig. [T5b. Capillary forces being attractive, 
is a tensile stress. For P* = 0.436, in the quasistatic 
limit, this contribution increases up to 20% in magni¬ 
tude. Consequently, the positive contribution of contact 
forces to reaches about 1.2cr^^ for P* = 0.436. 

The relative importance of the contributions of con¬ 
tacts and distant capillary forces to and is simi¬ 
lar; in the quasistatic limit and for P* = 0.436, one has 
~ —0.16 and ~ —0.25. 

Accordingly, the contribution of normal contact forces 
is the dominant one in normal stress differences A^i, N 2 
(with a notable contribution of tangential forces to N 2 , 
typically 20% at low P*). 

2. Elastic-frictional forces and capillary forces 

An alternative decomposition of the stress tensor is: 

CTa/3 — (14) 

in which is the contribution of capillary forces (either 
in the contacts or for distant interacting pairs), is 
the contribution of normal elastic forces and is the 
contribution of tangential forces. 

The normal elastic forces contribute more than 90% of 
the shear stress, whatever P* and I. 

The contribution of tangential forces to the normal (di¬ 
agonal) elements of the stress tensor is negligible, but 
that of capillary forces is very important: for P* = 0.436, 
negative terms (1 < o; < 3) are very large in mag¬ 
nitude: one observes < — 2cr„a for small P*, as 
shown in Fig. [T6| This large negative contribution is com¬ 
pensated by that of the repulsive normal elastic forces, 
^aa > 3(7aa- Such a large negative contribution of cap¬ 
illary forces to pressure implies that the particles are 


(a) 




FIG. 15: (Color online) Contribution of distant interactions 
to shear stress (a) and to normal stress (T22 (b). 

strongly pushed against one another, which increases the 
sliding threshold for tangential contact forces. 

Fig. [T71 shows the contribution of tangential forces to 
the total shear stress. As P* is decreased to P* = 0.436, 
the ratio increases to 0.18. Capillary forces con¬ 

tribute to the shear stress with the opposite sign (ct^p is 
positive while is negative). Fig. [171 shows that ratio 
cr“^/cri2 is always negative and decreases down to -0.12 
for P* = 0.436. Similarly to the case of normal stresses, 
the largest contribution is that of elastic normal forces, 
the (negative) contribution of which to compensates 
the (positive) term cr“P . 

D. Discussion 

One important clue to understand the enhanced shear 
strength of the cohesive material, as compared to the 
cohesionless, dry granular assembly, is the large tensile 
contribution of capillary force to normal stress: 

= -/ 3 ^ 


22 ’ 


(15) 
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FIG. 16: (Color online) Contribution of capillary forces to 
stress (722. 


with a coefficient j3 ranging, in the quasistatic limit, from 
about 0.15 (P* = 10) to 2.1 (P* = 0.436). Upon includ¬ 
ing the result for P* = 0.1 and / = 0.01, /3 reaches about 
7.2. This coefficient, and its variations with P*, can be 
approximately predicted from the values of solid fraction 
and coordination numbers. Contacts {zc^ on average, per 
grain) carry capillary force —Pq, and distant forces [zd 
per grain) average to a fraction of — Pq. Relation (IT^ can 
be used to evaluate the capillary contribution to pressure 

^zFci 

P, as- — < <- 7^- (This relation be- 

7ra^ 7ra^ 

tween P'=®'P and contact tensile strength Pq is sometimes 
referred to as the Rumpf formula, especiall y in the con¬ 
text of a prediction of rupture conditions |l3l . 

Dividing by , one obtains: 


P'=®'P <^Zc 

- < - < - 

ttP* 0-22 ttP* 


(16) 


Ignoring the small difference between P and , m pro¬ 
vides an estimate of coefficient (3 defined in (fT^ . Thus 
the value of (3 for reduced pressure P* = 0.436 is pre¬ 
dicted between 1.9 and 2.3 (and for P* = 0.1, it should 
reach about 8). Thus, quite unsurprisingly, the (nega¬ 
tive) relative contribution of capillary forces to normal 
stress if of order (1/P*) oc Pq/P, with a coefficient that 
may be deduced from $ and coordination numbers, ac¬ 
cording to CH). 

It is tempting to invoke a classical concept in geome¬ 
chanics, that of effective pressure, to describe the effect 
of capillary forces on the shear resistance of the material: 
the attractive forces create larger repulsive elastic reac¬ 
tions in the contact, corresponding to an effective pres¬ 
sure equal to (1 -I- /3)P. Furthermore, the local Coulomb 
condition in the contacts is to be written with those en¬ 
hanced normal repulsive forces. Capillary forces also con¬ 
tribute to shear stress, but, as apparent in Fig. 1171 in 
comparison to their influence on normal stresses, this is 
a small effect, and one may ignore it in a first approach. 



(b) 



FIG. 17: (Color online) Contributions of tangential (a) and 
capillary (b) interactions to total shear stress crj2. 


One assumes then that the shear behavior of the material 
is identical to that of a dry material under such effective 
normal stress cr®®. This approach leads to the follow¬ 
ing prediction for the P*-dependent quasistatic friction 
coefficient Pq: 


p*o = {l + f3)p^, (17) 


in which denotes the quasistatic internal friction coef¬ 
ficient for dry grains, P* = oo. Remarkably, if we further 
assume, as suggested by ([TC|). that [3 is roughly propor¬ 
tional to 1/P*, (3 ~ b/P*, we obtain a Mohr-Coulomb 
relation, Eq. [U for the stresses in the critical state: with 
the same value of internal friction as in the dry case, 
pI = p'iff and a macroscopic cohesion given by 


bp^Fo 


(18) 


Fig. [18] is a plot of 0-^2 versus CT22 ^ fbe yield locus - in 
which the predictions of relation [T3 both with the mea¬ 
sured coefficient f3 (Fig. [T51) . and with the one predicted 
as {z + Zc)^/{‘27tP*) from (flBl) . are confronted with the 
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FIG. 18: (Color online) versus (Tjj, in quasistatic flow, in 
units of Fo jo?. Square dots: numerical results (error bars are 
smaller); red crosses: predictions of (1171) . with exact coeffi¬ 
cient P\ blue circles: same with estimated /3. Insert: detail 
of numerical data for small P*, including additional point at 
P* = 0.1. 


numerical results. The admittedly crude prediction of 
relation (II3 appears surprisingly close to the numerical 
results on this plot. The relative error in the predic¬ 
tion for stress ratio with the measured value of /3, 
is actually about 5% at P* = 10, increasing to 20% at 
P* = 0.436, and the value of /Xg for P* = 0.1 ( ~ 1.6) 
from the measurements at I = 0.01 is largely overesti¬ 
mated, at 2.7. 

One may directly test for the validity of a Mohr- 
Coulomb relation to the data by fitting a linear form for 
the data of FigfTSl Given the error bars (which are small 
and do not appear on the graph), an attempted straight 
line fit through all 5 data points with P* > 0.43 in Fig.fTSl 
is unambiguously rejected by the standard likelihood cri¬ 
terion. A linear fit is (barely) acceptable upon ignoring 
the value P* = 0.436, yielding nl = 0.340 ± 0.001 and 
o^cjFo = 0.267 ± 0.005 for the Mohr-Coulomb parame¬ 
ters. From m the predicted apparent macroscopic co¬ 
hesion is above O.iFo/a^, and varies according to which 
data are used to identify b in (ITOl) . The result ~ 1.6 
for P* = 0.1 (corresponding to = fi*P* = 0.16) 

is thus in contradiction with the Mohr-Coulomb model, 
which becomes increasingly inadequate for smaller P*, 
as apparent in the insert in Fig. 1181 

The performance of the simple effective pressure pre¬ 
diction for the P* dependence of /ig is better visualized 
in Fig. [ini which, unlike Fig. [THl is not sensitive to stress 
scale. The global increase of /ig is predicted, yet overes¬ 
timated for the smallest P* values. 

One aspect that is not captured by this approach is 
the dependence of /ig on meniscus volume (Tab. lIVI) : the 
variations of coefficient /3 (from 1.8 to about 2.2 as V 
increases from 10“®a^ to 5.10“^a^) are insufficient to ac¬ 
count for the increase of the friction coefficient. 

There are quite a few reasons for the effective pres¬ 
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FIG. 19: (Golor online) Apparent quasistatic friction coeffi¬ 
cient /xq versus 1/P* - showing the value of /xo° for 1/P* = 0. 
Measurements and predictions of (HU. with exact and esti¬ 
mated coefficient /3, same symbols as in Fig. 1181 


sure approach to fail: while the mechanical properties 
are supposed to be the same once stresses are corrected, 
the density of the material, for one thing, is different in 
the dry and the wet case (with <i> varying between 0.525 
and 0.595 as P* grows from 0.46 to infinity); capillary 
forces also contribute to shear stress, the force network 
is bound to be different, etc. Nevertheless, although ad¬ 
mittedly crude, the prediction based on dn} proves apt 
to capture the trend of the change of /Xg with P*, al¬ 
though it overestimates its growth at small P*. As to the 
Mohr-Coulomb representation of yield stresses, it might 
be used as an approximation for P* > 1, but the obser¬ 
vations clearly preclude the definition of unique values of 
macroscopic cohesion and friction coefficient according 
to ([1]) for smaller pressures. 

Ref. reports on a laboratory study of quasistatic 
yield loci versus at the onset of plastic yielding 
and flow) of various kinds of wet granular assemblies in 
the pendular regime, including glass beads, which offer a 
suitable experimental comparison to our results. It pro¬ 
poses (under the name ’shift theory’), exactly the same 
effective pressure approach as the one we have attempted 
here, and concludes that it provides a good approxima¬ 
tion, by which the yield condition of wet materials is de¬ 
duced from the one of the dry grains. Interestingly, the 
investigated P* values in this study range from about 0.2 
to ^ 2.5, and the yield locus is slightly concave, as in our 
numerical results. Measured values of /x* are similar to 
our results (with e.g., /x* ~ 0.7 for P* = 1), and little 
change is obtained by increasing saturation by a factor of 
3. Some possible differences between those experiments 
and our simulations could result from the different state 
of the material (the experiments are not necessarily car¬ 
ried out in steady state quasistatic shear flow, and could 
depend on the initial assembling process), and from the 
value of the wetting angle. However, quite a satisfactory 
semi-quantitative agreement is obtained. 








15 


In the following sections, for a better assessment of 
the rheophysical effect of attractive capillary forces, mi¬ 
croscopic and microstructural aspects of force networks 
are investigated in greater detail. 


VI. FORCE DISTRIBUTION 

The distribution of intergranular force values in a gran¬ 
ular material in equilibrium [sil - ls^l . or in inertial 
flow 0,11 has received a lot of attention in the recent 
literature. While the probability distribution function 
of force values in cohesionless systems tends to decrease 
exponentially, on a scale given by the average (Tat), in 
cohesive granular assemblies, characterized by the con¬ 
tact tensile strength Fq, the equilibrium force distribu¬ 
tion evolves, as P* decreases to low values, towards a 
roughly symmetric distribution about zero, with values of 
both signs of order Fq (F/v) Blii. As compared to 
the two-dimensional results of Refs. liilil, the present 
3D numerical study of wet spherical grain assemblies does 
not investigate very small P* states, but involves longer- 
ranged distant interactions. The positive force wing of 
the p.d.f. of normal forces near the quasistatic limit is 
shown in Fig. 1201 showing the gradual departure from 



FIG. 20: (Color online) Distributions of normal forces normal¬ 
ized by average normal force (F^) for I = 10“® and different 
values of P*. 



FIG. 21: (Color online) Distributions of normal forces for 
P* = 0.436 and different values of I, normalized with Fq. 


results of Fig. [T^l a depletion of the population of con¬ 
tacts, compensated by a greater number of distant grains 
joined by a meniscus. To understand better the distribu¬ 
tion shape for negative values. Fig. [22] distinguishes the 
distributions of contact and distant (attractive) forces. 
Contact force distributions exhibit a maximum in zero, 
with negative values becoming more frequent as I in¬ 
creases. The larger value of the pdf near —Fq signals 
then the opening of more contacts. The distant interac- 



F^/Fq 


the cohesionless distribution shape, and the transition to 
a cohesion-dominated force network with values of order 
Fq, ratio (F/v)/Fo being approximately proportional to 
P* as discussed in Sec. IV Bl 

At low reduced pressure, as for P* = 0.436, it is more 
appropriate to normalize the distribution by Fq, as in 
Fig. 1211 This plot shows the influence of inertia param¬ 
eter I, which is, for large positive values, qualitatively 
similar to the one observed with dry grains; the distribu¬ 
tion widens, large forces being associated with collisions 
between grains or groups of grains. Another effect of 
increasing the inertial number is, as expected from the 


FIG. 22: (Color online) Gontributions of contact and distant 
interactions to p.d.f. of normal forces for P* = 0.436, and 
two values of /, 7 = 0.316 and I = 0.001. The vertical dashed 
line corresponds to force at rupture distance, F'^^^(F)o). 

tions are responsible for the non-monotonic part of the 
pdf. On the one hand, the sharp maximum near —Fq sig¬ 
nals a large population of grain pairs at close distance, in 
agreement with the fast increase of z(h) at small h visi¬ 
ble in Fig. [TH On the other hand, the increase near the 
minimum attractive force at rupture distance Dq merely 
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reflects the slow variation of function (Fig. [2]). 

The ” effective pressure” concept relies on the assump¬ 
tion that the effect of attractive capillary forces are simi¬ 
lar to that of a larger applied isotropic stress. One way to 
test such an idea at the microscopic scale is to compare 
the distributions of normal elastic forces: if normalized 
by the average elastic force, related to the effective pres¬ 
sure, those should be independent on P* and similar to 
the force distribution in a cohesionless system. Fig. [53] 



FNe/<FNe> 


FIG. 23: (Color online) Distributions of normal elastic forces 
at small I for different P* values. 

compares the distributions of elastic normal forces, nor¬ 
malized by their mean value, for small / and different 
values of P*. Those distributions are roughly similar, 
but show, as expected, notable discrepancies for values 
of order Fq. The decay for large values is faster in cohe¬ 
sive systems, reflecting a difference in force networks. 


VII. AGGLOMERATION 

The aggregation of cohesive grains is observed and re¬ 
ported in many numerical and experimental studies, and 
is exploited in industrial processes [H, It was di¬ 
rectly observed in flow of cohesive granular assemblies, 
both in numerical model materials [l^ . and in exper¬ 
iments with wet powders [5611. A numerical study of 
steady state chute flow reports an increase of the 
number of long-lasting contacts in the presence of cohe¬ 
sive forces. Weber et al. in carried out a detailed 
study of the effect of capillary forces on agglomerate du¬ 
ration and size. The agglomeration phenomena in steady 
shear flow is studied here, first, by measuring contact 
ages and meniscus ages, depending on state parameters. 
Then, the age-dependent size of clusters is measured, de¬ 
pending on P* and I. These clustering properties are 
related to the material rheology. 



FIG. 24: (Golor online) Distribution of the age of contacts 
for different values of P* and I — 0.1. Inset shows the same 
graph in a shorter range of r'^'y. 


A. Age of contacts and of distant interactions 

The distribution of the age of contacts for I = 10“^ 
and different values of P* is shown in Fig. [531 P{t'^'^) is 
the probability distribution of contact ages Tc, expressed 
as a strain The decrease of P(t'^ 7 ) is slower for 

smaller P *, showing that for the stronger cohesive forces 
the contacts survive over larger strain intervals . 

For large enough strains, > O.b, these probabil¬ 
ity density functions decay with an exponential form, 
P(r'^ 7 ) oc Values of decay times Tq, given in 

Tab. EH increase as we decrease P*. Average contact 
ages, also provided in the table, show the same be¬ 
havior (r^vg smaller than Tq because the distribution 
is not exponential for short times - see inset on the fig¬ 
ure). Fig. [55] shows the evolution of the pdf with / for 
two different values of P*, revealing, as expected, that 
contact ages (in units of I/ 7 ) decrease in faster flows. 
For I < 10“^, curves appear to coincide, showing nearly 
quasistatic behavior. The probability distribution func¬ 
tion of the age of interactions P{Pj) (i.e., the age of 
liquid bridges) is also shown in Fig. [55] for different val¬ 
ues of P*. Liquid bridges survive for quite large strain 
intervals, reaching several units with a probability of or¬ 
der 0.1, which increase as P* decreases. Initially, most 
liquid bridges survive at least for strains of order 0.1. Be¬ 
yond unit strain curves might be fitted by an exponential 
function too, defining a decay time Tq. P*-dependent 
values of Tq and of the average meniscus age are 
listed in Tab. ED Remarkably, the curves do not present 
any notable difference for different values of I: the pairs 
may lose their contacts in faster flows, but they are still 
bonded with liquid bridges. The age of contacts and of 
distant interactions thus reveal the formation of aggre¬ 
gates in the presence of capillary forces. These clusters 
are transported by the flow for some distance before they 
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FIG. 25: (Color online) Distribution of the age of contacts for 
different values of I and two different values of P*. 



FIG. 26: (Color online) Distribution of the age of menisci for 
different values of P* (same for all values of I). Inset shows 
the same graph in a shorter range of r‘7. 

are broken or restructured. They may survive for strain 
intervals of a few units. 

B. Clusters 

Clusters are defined as sets of grains connected by 
liquid bonds for a minimum time, and the cluster¬ 
ing tendency might be appreciated on recording the 
dependent mass-averages, cluster size: 

(19) 

in which the summations run over all clusters i containing 
Si grains. The results for for the different values 

of P* and I are shown in Fig. [271 as functions of 
For small values of almost all particles are gathered 


p* 

■yrS 

T 7avg 
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T 7avg 

0.436 

0.306 

0.258 

1.704 

1.609 

1 

0.180 

0.154 

1.437 

1.325 

2 

0.153 

0.111 

1.295 

1.187 

5 

0.128 

0.087 

1.164 

1.072 

10 

0.120 

0.080 

1.102 

1.021 

00 

0.118 

0.074 

— 
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TABLE VI: (Color online) Decay time of age distribution 
function for contacts, t§, and for all interactions, tq, obtained 
by an exponential fit to the data of Fig. 1241 and of Fig. I26II : 
average contact age r^vg and interaction age r^vg, for different 
values of P* and I = 0.1. All four times are normalized by 
shearing time I/7. 

in a single cluster, which results in {S^^)ra — 4000. For 
large values of r ‘^'7 most particles are isolated or part of 
a very small cluster, and so a value of order 1 for 
is expected. In the midrange of increase for 

decreasing P* or /, i.e. for stronger capillary forces or 
slower flows. 



FIG. 27: (Color online) Mass-averaged cluster size, 
versus cluster age for different values of P* and I. 


VIII. FABRIC ANISOTROPY 

The capacity of granular assemblies to form anisotropic 
force networks is the only origin of shear strength with 
frictionless grains [a j and is known to play a cen¬ 

tral role in the shear strength of frictional grains as well. 
To understand how contact and distant capillary forces 
contribute to the shear stress, it is instructive to study 
the distribution of contact orientations (normal vectors 
n on the unit sphere S), E{n), which, to lowest order, is 
characterized by the fabric tensor: 

Pa/3 — {‘^a'^/n') — / P(^Tl')n(yTl//j d Tl 
Jn 


( 20 ) 
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FIG. 28: (Color online) Sketch of approaching {ip near Stt/I) 
and receding pairs {ip near 7r/4) in macroscopic shear flow. 


The connection between fabric and normal force contri¬ 
bution to stresses, or , is quite direct, as one has: 

j E{n){F^)nnan/3(fn, ( 21 ) 

an integral over the unit sphere in which {F^)n denotes 
the average normal force carried by the pairs with ori¬ 
entation n. As reported in Sec. IV Cl the contribution of 
the normal forces, amounts to more than 80% of the 
shear stress. The contribution of fabric parameters F^^ 
to shear stress might be visualized in Fig. [28l On av¬ 
erage, if pairs are preferentially oriented with the normal 
vector within a compression quadrant in the shear flow, 
then Fi 2 < 0 will tend to increase the absolute value of 
if forces are positive, and to decrease it they are neg¬ 
ative. On the other hand, negative forces will increase 
the absolute value of if preferentially oriented in the 
extension quadrants in the shear flow. 

The evolutions of fabric parameters F^^ and F ^^, per¬ 
taining respectively to contact and distant normal forces, 
versus /, for different P* values, are displayed in Fig. 

F^^ is negative, signaling the contribution of normal con¬ 
tact forces to shear strength (as {F^'"^) is the dominant 
contribution to (F^), related to V by (fT^ ). The largest 
value, and the largest variation of F^^ with I, is ob¬ 
tained in the dry system [P* = oo). The decrease of this 
anisotropy parameter for smaller P* may be understood 
in reference to the clustering phenomena and to the larger 
duration of contacts evidenced in Sec. IVIII Longer-lived 
contacts rotate in the shear flow, and are less favorably 
oriented in the compression quadrant. Shear flow carries 
agglomerates for some distance before they break and 
thus their random tumbling motion increases the isotropy 
of the contact orientations. In faster flows (for larger 
/), while contacts tend to open in cohesionless materials, 
enhancing fabric anisotropy, cohesive contacts can resist 
flow agitation and inertial effects better, whence a smaller 
I influence; if they open, they transform into attractive 
distant interactions, and the anisotropy of distant inter¬ 
actions also decreases. Those distant capillary forces are 
characterized by a comparatively large anisotropy, about 
three times as large as |. As is positive, those dis¬ 
tant attractive forces contribute to increase the internal 
friction coefficient. Unlike F^^ increases for smaller 



I 

(b) 



I 

FIG. 29: (Golor online) Fabric parameter F 12 for contacting 
pairs (a) and distant interactions (b) versus I for different P*. 


P* values, which corresponds to the growing contribu¬ 
tion of distant interactions to shear strength shown in 
Fig-HSI The different rule of meniscus formation (at con¬ 
tact) and breakage (at distance Dq) explains, in part, this 
large fabric anisotropy of attractive forces: approaching 
particles are not attracted to each other, whence a small 
number of distant interacting pairs in the compressive 
quadrant, with a negative contribution to F 12 ; as par¬ 
ticles get separated, receding pairs are still attracted to 
each other, whence a positive contribution to F 12 from 
the extension quadrant. In the model without menis¬ 
cus hysteresis, assuming capillary attraction to appear as 
soon as grains approach within distance Dq, strongly 
decreases, from 0.14 to about 0.07 at P* = 0.436 and 
small I. 


IX. SUMMARY AND DISCUSSION 

The rheological properties of unsaturated granular ma¬ 
terials, in which a small amount of wetting liquid, form- 
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ing liquid bridges and transmitting attractive capillary 
forces between particles, generalize, in many respects, 
previous observations on cohesive granular materials, 
with macroscopic properties exhibiting similar depen¬ 
dences on I and P*. Thus, compared to dry materials, 
the apparent internal friction coefficient = aril 012 
is enhanced (from 0.33 to more than 1 in the explored 
range P* > 0.1); looser structures are stabilized, even 
in the quasistatic limit ($ ~ 0.52, for P* = 0.436 is be¬ 
low all packing densities with cohesionless grains), even 
though contact coordination numbers, due to the absence 
of rattlers, may be larger. Our results describe those ef¬ 
fects in quantitative form, in the range P* >0.1 and 
10“^ < I < 0.56, and specify the dependence on various 
features of the model. We only predict quite a small de¬ 
pendence of the rheology on saturation within the pendu¬ 
lar range (up to 5-10%), in agreement with experimental 
observations [l^[50l|. 

More accurate models of the capillary force dependence 
on intergranular distance, or of the distribution of liq¬ 
uid between menisci with varying volumes, would hardly 
change the results. Interestingly, though, some variants 
of the model, although arguably not realistic, have no¬ 
tably different rheological properties. Thus, reducing the 
meniscus volume to very small values would have quite 
a notable effect on internal friction and density - but, 
in practice, menisci are unlikely to form with such small 
liquid contents. Assuming menisci form as soon as grains 
approach to their maximum extension distance (range of 
capillary force) would also strongly affect macroscopic 
properties. 

Shear localization systematically affects shear flows at 
low P*, and we could not measure the constitutive be¬ 
havior at P* =0.1 except for some intermediate / val¬ 
ues of order 0.01. Localized states are characterized by 
velocity profiles with gradients concentrated within nar¬ 
row bands, where the solid fraction is well below its bulk 
value. The band thickness lies in the range of 5 to 10 
grain diameters a at small I, but might be as small as 
about 1.5a in faster flows (I > 0.1). 

We also record normal stress differences, which are 
larger than for dry grains, and tend to grow with de¬ 
creasing P*. The second normal stress difference, in par¬ 
ticular, reaches 20% of the imposed normal stress a 22 in 
the quasistatic limit for small P*. 

The effective pressure approach to the yield criterion 
of wet grains ignores such sophistications, as well as den¬ 
sity or microstructural changes due to capillary forces. 
It assumes critical states to be in correspondence for dif¬ 
ferent values of P*, as though the introduction of capil¬ 
lary forces, pushing grains against their neighbors, were 
equivalent to the application of a larger confining pres¬ 
sure. Such a crude approach is in fact surprising suc¬ 
cessful as a rough approximation, to predict the increase 
of fi* for decreasing, but not too small, P*, say P* > 1 
(below P* = 1, the increase of /r* is overestimated). The 


effective pressure might be evaluated upon adding, to the 
applied pressure, the capillary contribution to the aver¬ 
age stress. This contribution might itself be estimated 
from density and coordination numbers, which leads to 
a Mohr-Coulomb form (HD for the variation of tensile 
strength p^(J 22 with normal stress a 22 (if 4* and z do not 
vary too much). Such a form of the critical state plas¬ 
ticity criterion proves however inadequate to describe the 
whole range of reduced pressures P*: data are incompat¬ 
ible with a P*-independent macroscopic cohesion c. 

Another remarkable feature of the measured rheology, 
compared to dry granular materials, is the slow variation 
with I of both macroscopic rheological parameters such 
as /X*, and microscopic data such as coordination num¬ 
bers. Although we could not observe a direct correlation, 
a slower increase of function /i*(/) could signal a shear 
banding tendency. 

Many of those rheophysical features are explained by, 
or, at least, related to, the strong clustering tendency 
emerging as attractive forces gradually become the dom¬ 
inant ones, upon decreasing P*. As contacts are stabi¬ 
lized by attractive forces, they do not so easily open as 
the network is being sheared. When they do (which hap¬ 
pens preferentially in the extension direction within the 
average shear flow, whence a fabric anisotropy of distant 
interactions contributing to shear strength), the network 
of grains bonded by liquid bridges might still be con¬ 
nected, forming enduring connected clusters. The sur¬ 
vival of such clusters over quite notable strain intervals 
(reaching several unities with sizable probability) should 
limit the dilating tendency of faster flows. It also main¬ 
tains a network in which the capillary forces act in closer 
similarity to an effective pressure. 

Admittedly, this is still a descriptive rheophysical sce¬ 
nario. More quantitative studies should be carried out 
to better characterize the deformation mechanisms of 
the grain clusters. The shear banding phenomenon cer¬ 
tainly deserves detailed investigations, in which sample 
size and shape effects should be systematically assessed, 
and partly localized velocity an density fields analyzed 
and related to a stability analysis. 

Although we argued that our measurements agree with 
the limited available experimental results, more labora¬ 
tory data should certainly be used, with, if possible both 
information on rheology and on micromorphology and 
liquid distribution. 

Finally, a complete numerical model should eventually 
be designed, capable of dealing with saturations exceed¬ 
ing the limited pendular range, and of describing the liq¬ 
uid motion. One should thus be equipped to model the 
mixing process as well as the rheology of the homoge¬ 
neous mixture of the grains and the liquid. The Lat¬ 
tice Boltzmann method for a diphasic interstitial fluid 
medium, coupled to a DEM description of grain motion, 
is a promising perspective [13, [6l|. 
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